/*
 * 蒙特卡洛重要抽样积分： \int^50_0 (3-x)*exp(-x) dx
 */

#include<iostream>
using namespace std;
#include<cmath>

#include"library.h"

double f(double x){
	return (3-0.1*x)*exp(-x);
}

double g(double x){
	return exp(-x);
}

/*
 * 蒙特卡洛积分：简单抽样
 */
void simple_sampling(double xa, double xb, double (*f)(double x), int n){

	int i;
	double x, y, sum_y = 0, sum_y2 = 0, ave_y, ave_y2, I1, I2;

	for(int i=0;i<n;i++){
		x = (double)rand()/RAND_MAX*(xb-xa) + xa;//取随机数 x 在 [xa,xb] 均匀分布
		y = f(x);
		sum_y += y;
		sum_y2 += y*y;
	}
	ave_y = sum_y / n;
	ave_y2 = sum_y2 / n;

	I1 = (xb - xa) * ( ave_y - 4 * sqrt( (ave_y2 - ave_y * ave_y)/n ) );//置信区间左边界
	I2 = (xb - xa) * ( ave_y + 4 * sqrt( (ave_y2 - ave_y * ave_y)/n ) );//置信区间右边界
	cout<<"Simple sampling: integration value is "<< (I2+I1)/2<<" +- "<<(I2-I1)/2<<" by "<<99.993666<<"% chance."<<endl;
	// 4 sigma 置信度为 99.993666%
}

/*
 * 蒙特卡洛积分：重要抽样，抽样概率密度为 g(x)
 */
void importance_sampling(double xa, double xb, double (*f)(double x), double (*g)(double x), double gmax, int n){

	int i;
	double x, y, sum_y = 0, sum_y2 = 0, ave_y, ave_y2, I1, I2;

	for(int i=0;i<n;i++){
		x = gnr_rand(xa, xb, g, gmax);
		y = f(x)/g(x);
		sum_y += y;
		sum_y2 += y*y;
	}
	ave_y = sum_y / n;
	ave_y2 = sum_y2 / n;

	I1 =   ave_y - 4 * sqrt( (ave_y2 - ave_y * ave_y)/n ) ;//置信区间左边界
	I2 =  ave_y + 4 * sqrt( (ave_y2 - ave_y * ave_y)/n ) ;//置信区间右边界
	cout<<"Importance sampling: integration value is "<< (I2+I1)/2<<" +- "<<(I2-I1)/2<<" by "<<99.993666<<"% chance."<<endl;
	// 4 sigma 置信度为 99.993666%
}

int main(){

	int n=1E6;
	simple_sampling(0, 50, f, n);
	importance_sampling(0, 50, f, g, 1, n);
	cout<<"exact value = 2.9"<<endl;

	return 0;
}
